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^ Abstract 

We study the special point in the phase diagram of a semi-infinite system, where the bulk 
transition is in the three-dimensional Ising universality class. To this end we perform a finite size 
scaling study of the improved Blume-Capel model on the simple cubic lattice with two different 
types of surface interactions. In order to check for the effect of leading bulk corrections we have 
also simulated the spin- 1/2 Ising model on the simple cubic lattice. We have accurately estimated 
the surface enhancement coupling at the special point of these models. We find yt^ = 0.718(2) and 
yfii = 1.6465(6) for the surface renormalization group exponents of the special transition. These 
results are compared with previous ones obtained by using field theoretic methods and Monte 
Carlo simulations of the spin-1/2 Ising model. Furthermore we study the behavior of the surface 
transition near the special point and finally we discuss films with special boundary conditions at 
one surface and fixed ones at the other. 
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I. INTRODUCTION 



In this paper we shall study the special point in the phase diagram of a semi-infinite 
system. For reviews on surface critical phenomena see refs. [ll-j^. Let us briefiy recall the 
basic features of this phase diagram at the example of the spin- 1/2 Ising model on the simple 
cubic lattice and a semi-infinite geometry. It's reduced Hamiltonian is given by 

H = -(3 s^Sy- h'^s^- I3i ^ s^Sy-hi'^s^ , (1) 

<xy> X <xy>^S xGS 

where x and y denote sites of the lattice, < > is a pair of nearest neighbors and S is 
the surface of the system. The spin Sx at the site x can take either the value —1 or 1. 
The Boltzmann factor is given by exp(— if), since the reduced Hamiltonian incorporates the 
temperature. We define (3 = J/ksT, (3i = Ji/ksT, where J is the coupling constant in the 
bulk, Ji the excess coupling constant at the surface and T is the temperature. Below we 
shall refer to (3 and f3i as the coupling in the bulk and the excess coupling at the surface, 
respectively. Below we shall consider vanishing external fields h = hi = in the bulk as well 
as at the surface. 

In figured] we have sketched the phase diagram of this system. For /3 > /3c, where (3c is the 
critical coupling of the bulk system, the spins in the bulk are ordered. As a consequence, also 
the spins at the surface are ordered. At /3 = the spins at the surface decouple completely 
from those of the bulk. Hence a two dimensional Ising model remains that undergoes a 
phase transition at f3i = /3c,2D- Starting from the point {0,(3c,2d) there is the line of surface 
transitions, where the spins at the surface order, while those of the bulk remain disordered. 
This line hits the line [Pc, Pi) in the so called special or surface-bulk point, which is a tri- 
critical point. The transitions from disordered surface and disordered bulk to ordered bulk 
and ordered surface are called ordinary transitions, while those from disordered bulk and 
ordered surface to ordered bulk and ordered surface are called extraordinary transitions. 

Surface critical phenomena had been studied first by using the mean-field approximation 
[l|. The application of field theoretic methods to surface critical phenomena is complicated 
by the fact that translational invariance is broken by the surface. As a consequence, surface 
critical exponents are only computed up to O(e^) in the e-expansion j^, [3|. Furthermore 
surface critical phenomena have been studied by using high temperature series expansions, 
real space renormalization group methods and Monte Carlo simulations of lattice models. 

The special point is characterized by the two relevant bulk renormalization group (RG)- 
exponents yt and yh and the two relevant surface RG-exponents yt^ and yt^- Similar to the 
pure bulk case, these RG-exponents can be related to a number of surface critical exponents 
that characterize the behavior of thermodynamic quantities related to the surface in the 
neighborhood of the special point [l|-[3|. 

In the present work we locate the special point of three different lattice models and 
determine the RG-exponents yt^ and yhi by using finite size scaling (FSS) Q methods. In the 
presence of a surface, typically corrections oc appear [iMsl], where L is the linear extent 
of the finite system. Analysing numerical data, it is difficult to disentangle these corrections 
from leading bulk corrections which are oc L-'^, where lo = 0.832(6) Q. Therefore, in 
addition to the Ising model we simulate the Blume-Capel model that is a generalization of 
the Ising model. In addition to —1 and 1 the spin can take the value 0. The parameter D 
of the Blume-Capel model controls the density of spins with Sx = 0. It has been shown that 
for a particular value D* of the parameter D the amplitude of leading corrections to scaling 
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FIG. 1. Sketch of the phase diagram of the semi- infinite system. On the x-axis we plot the 
couphng P of the bulk and on the y-axis the excess coupling f3i of the surface. The phase, where 
both the boundary spins and those of the bulk are disordered is labelled by A, while the one with 
disordered bulk and and ordered surface is labelled by B. The phase, where the bulk and the surface 
are ordered is labelled by C. The line of surface transitions is labelled by S, the line of ordinary 
transitions by O and the line of extraordinary transitions by E. These three lines meet in the so 
called special or surface-bulk point that we have labelled by SP. A discussion is given in the text. 



vanishes. See [5| and refs. therein. The precise definition of the model is given below in 
section [Tll 

The outline of the paper is the following: First we define the models that we simulate. 
We summarize results for the critical coupling and the critical exponents of the bulk system. 
We define the quantities that we have measured and discuss their finite size scaling behavior. 
The scaling behavior is affected by corrections, which have to be taken into account when 
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analysing Monte Carlo data. Then we report our numerical results: First we have simulated 
three different models on Lr" lattices, where L is the linear extend of the lattice. Based 
on these simulations we determine Pi^sp for these models and obtain estimates of yt^ and 
Next we study the surface transition in the neighborhood of the special point. Then we 
discuss the magnetisation profile of films with special boundary conditions at one surface and 
fixed boundary conditions at the other. Finally we summarize and compare our results for 
the RG-exponents with those obtained by field theoretic methods and previous simulations 
of the spin-1/2 Ising model. 

II. MODELS AND OBSERVABLES 

The Blume-Capel model on the simple cubic lattice is characterized by the reduced Hamil- 
tonian 

<xy> X X 

where x = (xq, xi, X2) denotes a site of the lattice. The components xq, xi and X2 take integer 
values in the range 1 < Xj < Lj. The spin might take the values —1, or 1. In the following 
we shall consider a vanishing external field /i = throughout. The parameter D controls 
the density of vacancies = 0. In the limit D — )■ — 00, these vacancies are completely 
suppressed, and hence the spin-1/2 Ising model is recovered. For —00 < D < Dtri the 
model undergoes a second order phase transition in the three-dimensional Ising universality 
class. For D > Dtri the transition is of first order. The most recent estimate for the tri- 
critical point is Dtri = 2.0313(4) [6|. Numerically, using Monte Carlo simulations it has 
been shown that there is a point [D* , Pc{D*)) on the line of second order phase transitions, 
where the amplitude of leading corrections to scaling vanishes. Our most recent estimate is 
D* = 0.656(20) [5|. In [sl we have simulated the model ai D = 0.655 close to Pc on lattices 
of a linear size up to L = 360. From a standard finite size scaling analysis of renormalization 
group invariant quantities such as the Binder cumulant we found 

;ge(0.655) = 0.387721735(25) (3) 

for the critical coupling of the bulk at D = 0.655. Recent estimates for the critical coupling of 
the bulk of the spin-1/2 Ising model are /3c = 0.22165455(3), table X of [7] and 0.22165463(8) 
[sl. In the following we assume (3c = 0.2216546. The amplitude of leading corrections to 
scaling a.t D = 0.655 is at least by a factor of 30 smaller than for the spin-1/2 Ising model. 

Our recent estimates for bulk critical exponents in the three-dimensional Ising universality 
class are j5|] 

z/ = 0.63002(10) , (4) 
r/ = 0.03627(10) , (5) 
w = 0.832(6) . (6) 

A. Film geometry and boundary conditions 

In Monte Carlo simulations we are restricted to finite lattices. Therefore a surface requires 
a counterpart. This means that we actually study systems with a film geometry. In the ideal 
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case, film geometry means that the system has a finite thickness Lq, while in the other two 
directions the thermodynamic limit Li, L2 — )■ oo is taken. In order to approximate this limit 
in Monte Carlo simulations, one usually chooses Lq ^ Li, L2 and applies periodic boundary 
conditions in the 1 and 2 directions. Note that we shall simulate lattices with Li = L2 = L 
throughout. As we shall see below, in order to compute surface critical exponents, the 
condition Lq <^ Li,L2 is not mandatory. Our estimates for the surface critical exponents 
are actually obtained from simulations of lattices with Lq = Li = L2 = L. 
The reduced Hamiltonian of the Blume-Capel model with film geometry is 

H = -(3 s^Sy + D'^sI- h'^s^- I3i ^ s^Sy-(32 ^ s^Sy 

<xy> X X <xy>,xo=yo=l <xy>,xo=yo=Lo 

+ DiY,sl + D2 si - hi Sx - h2 ' 

x,xq=1 x,xo=Lo x,xo=1 x,xq=Lq 

where we have put the surfaces at xq = 1 and xq = Lq. In our convention < xy > runs over 
all pairs of nearest neighbor sites with fiuctuating spins. Note that here the sites (l,xi,a;2) 
and {Lo,xi,X2) are not nearest neighbors as it would be the case for periodic boundary 
conditions. For each bulk term there is a corresponding surface enhancement term. In our 
simulations and the analysis of the data we have used the surface couplings 

^i = f3i + f3 and ^2 = 1^2 + /3 (8) 

as parameters instead of the excess surface couplings /3i and (32- 

Note that as long as the bulk transition and the line of surface transitions remain continu- 
ous, the qualitative features of the phase diagram that we have discussed in the introduction 
should remain unchanged. Since the values of D, D + Di and D + D2 that we shall consider 
are much smaller than Dtri = 2.0313(4) [6| and Dtri,2D = 1.966(2) [sl] of the two-dimensional 
Blume-Capel model on the square lattice, this should be the case in our study. 



B. Observables 

Renormalization group invariant quantities are very useful to locate critical or multicrit- 
ical points. We study the Binder cumulant 

where m = Ylx ^x- The second moment correlation length is given by 



^2nd - \l — 77 5 (10) 
4sm tt/L 



where 

2ttxi • 



^explz— ^ ) (11) 



LoL2 \|^ -V Li 

\ X ^ 

is the Fourier transform of the correlation function at the lowest non-zero momentum in 1 
or 2 direction and x = {'^'^) / {.LqL"^) is the magnetic susceptibihty. Since in our simulations 
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Li = L2, the expectation value of F is identical for the 1 and the 2 direction. In order 
to reduce the statistical error, we have measured F for both directions and have averaged 
the results. The ratio C,2nd/L is renormalization group invariant. The third renormalization 
group invariant quantity that we consider is the ratio Za/Zp of partition functions, where Za 
is the partition function of a system with anti-periodic boundary conditions in 1 direction 
and periodic ones in 2 direction or vice versa, while in the case of Zp periodic boundary 
conditions are imposed in 1 and 2 direction. Also here, since Li = L2 in our simulations, 
we determine Za/Zp for both choices and average the results. The ratio Za/Zp of partition 
functions can be efficiently evaluated using the boundary flip algorithm |9|. Here we use a 
modified version of the boundary fiip algorithm as discussed in appendix A 2 of ref. jlo| . 
In the following we shall refer to the renormalization group invariant quantities U4, Za/Zp 
and C,2nd/L using the symbol R. For vanishing symmetry breaking fields h = hi = h2 = 0, 
Di = D2 and /3i = /32, a renormahzation group invariant quantity behaves as 

/?(Lo,L,/3,/3i) = g(L/Lo,t[Lo/eo]^%ti[^o/^i,o]^*0 , (12) 

where t = I3c — /3, ti = f3i^sp — f^i, ^0 is the amplitude of the correlation length in the high 
temperature phase and the normalisation factor ^10 is still undetermined. In section (lllip 
we have fixed the ratio L/Lq = 1 and have set the bulk coupling to its critical value: /3 = Pc- 
This allows us to determine the location of the special point /3i^sp by using the standard 
crossing method, where R is considered as function of The behavior of the slope of the 
renormalization group invariant quantities allows us to determine the surface RG-exponent 
yt^: Taking the derivative with respect to f3i we get 

dR{Lo y,Pi) ^ Q^L/Lo,t[Lo/^or,h[Lo/^i,or^) [^0/6,0]^^*^ (13) 
opi 

where Q is minus the partial derivative of q with respect to its third argument. 

Finally we define the surface susceptibilities for a vanishing surface magnetisation (mi): 

Xii = ^^ = ^'(r«?) (14) 

where 

rni = 

and 



Xl,X2 



Xi2 = ^ = L'{m,m2) (16) 



dh 



2 



where 



= 75 XI ^(■^0,a;i,X2) • (17) 



L2 

XI, X2 



The finite size scaling behavior of these quantities can be inferred from the singular part of 
the reduced free energy per area of the film. Its scaling form is 

/,(...,...,..., /ii, /i2) = Lo-"+^(7(Vi^o,. /ii[W4J^'^S/i2[V^hJ^'^2) , (18) 

where the constants //^ and remain undetermined here, i/hi and yh2 are the RG-exponents 
related to the surface fields. The susceptibilities defined above can be expressed as second 
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derivatives of fs with respect to the surface fields. Taking these derivatives on the right 
hand side of eq. (fT8|) we arrive at 

Xn = ^ = cL,'^^'[Lo/W^^ (19) 

and 

Xi2 = = cL^'^' [LoM''^ [L,W- . (20) 

In our study below, section [TlIBt y^^ = i/h^ and Ih^ = lh2, since we have chosen /3i = /32, 
Di = D2 and /ii = /12 = 0. 



C. Corrections to scaling 

Finite size scaling laws such as eqs. f lT2|13|19|l20|) are affected by corrections. These are 
caused by irrelevant scaling fields, the analytic background in the reduced free energy per 
area and the fact that e.g. in eq. f lT2|) the arguments of the scaling function q should be 
actually analytic functions of t and ti that we have linearized here. 

The leading bulk correction is oc Lq'^ , where our most recent estimate oj = 0.832(6) 
obtained from Monte Carlo simulations of the Blume-Capel and the Ising model jll is slightly 
larger than that obtained from field theoretical methods, e.g. u = 0.799(11) by using 
perturbation theory in three dimensions fixed and u = 0.814(14) by using the e-expansion 
There are also corrections oc Lq^, where n = 2, 3, ... . In the case of improved models 
these are highly suppressed and therefore we shall ignore them in the analysis of our data 
below. Following [l2| subleading corrections are characterized by C02 = 1.67(11). There is 
no reason to assume that the amplitude of the subleading correction vanishes in the case of 
the improved Blume-Capel model. Likely it is of similar size as in the spin-1/2 Ising model. 
Furthermore there are well established corrections with ~ 2, for example related to the 
breaking of the spatial rotational invariance by the simple cubic lattice |l3| . 

In the presence of surfaces there are also corrections caused by irrelevant surface fields. 
One expects that the leading ones are cx: Lq ^ jl|-[3[. 

In general we expect that for finite Lq finite size scaling laws such as eqs. (I13fl9f20p can 



be written in the form of a Wegner-expansion [14 



A(Lo) = cL5[l + aLo'^ + «*^o""1 (21) 

i=2 

with an infinite number of correction terms. Fitting data obtained from Monte Carlo sim- 
ulations, only a very limited number of terms can be taken into account. This unavoidable 
truncation of the Wegner-expansion leads to a systematic error of the estimate of e.g. the 
exponent x that is often larger than the statistical one. 

In the present work, we shall use ansaetze that include either no correction, a correction 
oc Lq ^ or corrections oc Lq ^ and oc Lq ^. In the case of the surface susceptibilities we shall 
also take into account a term for the analytic background. The systematic errors of our 
results are then estimated from the variation of the results obtained by using these different 
ansaetze. 
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TABLE I. Number of measurements divided by 10 . For each measurement we perform one heat- 
bath sweep, a certain number of single cluster and two wall-cluster updates. 



L 


BCl 


BC2 Ising 


8 


1000 


1000 


1000 


12 


1000 


1000 


1000 


16 


1000 


1000 


1000 


24 


1000 


1178 


877 


32 


1020 


1042 


718 


48 


861 


705 


305 


64 


593 


482 


224 


96 


301 


200 


132 


128 


202 


140 


116 



III. SIMULATIONS OF LATTICES AT THE SPECIAL POINT 

In this part of our study we consider three different models. In all cases we chose Lq = L, 
Pi = ^2-, Di = D2 and hi = h2 = 0. We have simulated the Blume-Capel model at 
D = 0.655 with two different choices of Di. The choice Z)i = is called BCl model and 
-Di — i- — 00 is called BC2 model in the following. Note that in the case of the BC2 model, 
the spins at the surfaces can take only the values —1 or 1. In addition, we have simulated 
the spin-1/2 Ising model. We have simulated at the best estimates of the bulk critical point; 
i.e. /3c = 0.387721735 for the Blume Capel model at D = 0.655 and (3c = 0.2216546 for 
the spin-1/2 Ising model. The numerical uncertainty of these numbers is negligible in the 
present study. 

For each measurement of the observables we performed the following sequence of Monte 
Carlo updates: First one sweep with a local update. In the case of the Ising model we use 
a local Metropolis and for the Blume-Capel model a local heat-bath update. One sweep 
means that we run through the lattice once in type-writer fashion. Then we performed a 
certain number of single-cluster updates [lH followed by two wall-cluster updates [l6j; one 
for each of the directions with periodic boundary conditions. In all our simulations we have 
used the Mersenne twister algorithm jl3] as pseudo-random number generator. 

In the case of the Blume-Capel models BCl and BC2, no previous estimate of the surface 
coupling /Si^sp was available. Therefore, we have successively improved our estimate of f3i^sp 
with increasing lattice size L. In order to obtain the observables as a function of (3i in 
the neighborhood of the simulation point, we have computed the coefficients of the Taylor- 
expansion of the observables in /3i around the simulation point up to third order. In tabledwe 
summarize the lattice sizes that we have simulated at and the statistics of these simulations. 
In total we have spent about 12 years of CPU time on a single core of a Quad-Core AMD 
Opteron(tm) Processor 2378 running at 2.4 GHz. 

In our simulations we have measured the renormalization group invariant quantities 
Za/Zp, Ui and i2nd/L. Analysing the data, it turned out that corrections to scaling are 
considerably larger for ^2nd/L than for Za/Zp and U^. Therefore we restrict the following 
discussion to Za/Zp and t/4. 
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TABLE 11. Fitting the ratio of partition functions Za/Zp for the BCl model with the ansaetze <I22|) . 
(HI or (ESP. 





/3l,sp,_BCl (^a/^p)* C 


d 




[22] 32 


0.5491558(23) 0.31318(5) 




14.2/3 


[22] 48 


0.5491482(32) 0.31341(8) 




2.32/2 


M 16 


0.5491430(33) 0.31374(10) 0.0131(14) 




6.08/4 


M 24 


0.5491350(47) 0.31406(17) 0.0194(31) 




0.69/3 


[55] 12 


0.5491423(49) 0.31375(19) 0.0121(47) 


0.023(34) 


8.87/4 


[25] 16 


0.5491280(69) 0.31445(32) 0.035(10) 


-0.216(94) 


0.51/3 



In a first step of tlie analysis we have computed tlie surface coupling /3i of the special 
point and the fixed point values of Za/Zp and f/4. To this end we have fitted the data with 
the ansatz 

R{^i,sp,L) = R* (22) 
where /9i,sp and R* are the free parameters of the fit and 

RiPi^sp, L) = R{l3,,s, L) + cSi,sp - hs) + j^ ihsp - hsf + j^ {hsp - hsf (23) 

where R{/3i^s,L), ci, C2 and C3 are obtained from the simulation at (3i^s- To check for the 
possible effect of corrections to scaling we have also used the ansaetze 

R{^i^sp,L) = R* + cL-' (24) 

and 

R{^i,sp, L) = R* + cL~^ + dL-^ . (25) 

First we have analyzed the three models separately. Let us first look at the results for the 
BCl model and the ratio Za/Zp. A selection of results is given in table [Til In these fits 
we have taken into account the data for all lattice sizes L > Lmin- Fits with an accept- 
able x^/DOF are obtained starting from Lmin = 48, 24 and 16 for the ansaetze f[22|) . f[M]) 
and f[2Sl) . respectively. Note that the differences of the results for (3i^sp and (Za/Zp)* for 
different ansaetze with an acceptable x^/DOF are larger than the statistical errors. Next 
we have fitted the Binder cumulant with the ansaetze f[221) . and f[2^ . In the case of 
the ansatz ([22]) we find xVDOF= 1.32/4 already for L^i^ = 24 with = 0.5491393(32) 
and Ul = 1.52338(7). Fitting the data for Lmin = 8 with the ansaetze ([24]) and ([25]) 
we get x^/DOF= 8.69/6 and 8.61/5, respectively. The results for the fit parameters are 
= 0.5491470(37), f/4 = 1.52303(10) and = 0.5491455(63), U4 = 1.52309(23), for 
the ansaetze and (125]) . respectively. Notice that the results for Pi^sp are compatible with 
those obtained from the analysis of Za/Zp. 

Next we have analyzed the data for the BC2 model. In the case of the ansatz fl22p larger 
Lmin are need than for the BCl model to get an acceptable x^/DOF. Fitting the data with 
the ansaetze fl24p and fl25p we see that the correction amplitude c is clearly larger for the 
BC2 than for the BCl model. 

Next we performed joint fits for the BCl and the BC2 model using the ansaetze ([Ml) and 
(125]) . In these fits we impose that, following universality, (Za/Zp)* and are the same for 



9 



TABLE III. Fitting the ratio of partition functions Za/Zp for the BCl and BC2 model with the 
ansaetze (f2i)) or ^2^. In the case of ansatz <r25)) we do not report d\ and d2 to keep the table 
readable. 

Ansatz L Pi,sp,BCi h,sp,BC2 {Zg/Zp)* ci C2 X^/DQF 

[21 12 0.5491418(20) 0.2940123(19) 0.31380(5) 0.0144(6) -0.0782(6) 19.78/11 
[H 16 0.5491461(24) 0.2940165(23) 0.31364(7) 0.0118(10) -0.0806(10) 9.36/9 
[21 24 0.5491416(35) 0.2940124(33) 0.31382(13) 0.0152(23) -0.0771(22) 6.20/7 
[25] 8 0.5491493(27) 0.2940202(26) 0.31346(9) 0.0045(17) -0.0863(16) 15.44/11 
[25] 12 0.5491481(28) 0.2940189(36) 0.31352(15) 0.0064(38) -0.0850(37) 15.13/9 
[25] 16 0.5491364(52) 0.2940086(49) 0.31407(23) 0.0242(70) -0.0664(68) 4.70/7 



TABLE IV. Same as table lllll but for the Binder cumulant instead of Za/Zp. 



Ansatz L h,sp,BCi Pi,sp,BC2 




16 0.5491396(42) 0.2940168(40) 1.52327(13) -0.005(2) 0.347(2) 16.75/9 

24 0.5491377(61) 0.2940104(57) 1.52345(23) 0.002(4) 0.348(4) 1.95/7 

8 0.5491520(45) 0.2940194(41) 1.52285(15) -0.011(3) 0.323(3) 13.58/11 

12 0.5491517(66) 0.2940202(62) 1.52284(27) -0.012(7) 0.324(7) 13.00/9 

16 0.5491342(82) 0.2940047(80) 1.52372(38) 0.017(12) 0.353(11) 5.84/7 



both models. A selection of results is given in table [XTTl We see that the correction amplitude 
C2 of the BC2 model is much larger than ci of the BCl model. We performed similar fits for 
the Binder cumulant f/4. Also here we observe that the correction oc has a much larger 
amplitude for the BC2 model than for the BCl model. The results for the surface couplings 
at the special point f3i^sp,Bci and Pi^sp,BC2 obtained from the analysis of and Za/Zp are 
compatible among each other. 
As our final result we quote 

{Za/ZpY = 0.3138(5) , Ul = 1.5234(10) (26) 

and 

^i,sp,Bci = 0.54914(2) , l3,,sp,BC2 = 0.29401(2) (27) 

where the error bars are chosen such that the results of all fits given in table IIIII and [IV] are 
covered. 

Finally we have fitted our data for the spin-1/2 Ising model using the ansatz (jUj). A 
selection of results is given in table |V] For Lmin = 16 and 24 we find acceptable values for 
X^/DOF. Nevertheless the result for (Za/Zp)* is not compatible with that obtained from 
the fits for the BCl and BC2 models. This could be explained by the fact that corrections 
oc L~'^ with Lj = 0.832(6) [51] are not explicitly taken into account. 

Finally we performed a number of different fits, where we take the results f[26|) for (Za/Zp)* 
and t/4 as input. In these fits we also include explicitly corrections cx L^'^. Taking into 
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TABLE V. Fitting the ratio of partition functions Za/Zp for the spin-1/2 Ising model with the 
ansatz ([HP. 



L hsp,i {Zg/ZpY c xVdof 

12 0.3330388(15) 0.31221(5) 0.00219(5) 9.03/5 
16 0.3330365(19) 0.31232(7) 0.00236(9) 4.34/4 
24 0.3330338(28) 0.31246(13) 0.00262(23) 2.77/3 



TABLE VI. Comparison of our result for the ratio Pi^spj/ Pc with those given in the literature. 



Ref. year Pi^sp,i/Pc 
[18] 1984 1.50(3) 
\19\ 1990 1.52(2) 
[20] 1993 1.5004(20) 
[21] 2005 1.50208(5) 

Here 2011 1.50243(9) 



account these different fits we arrive at the final estimate 

Pi,sp,i = 0.33302(2) (28) 

for the spin-1/2 Ising model. In table I VI I we compare our result with those given in the 
literature. Our result is compatible within error bars with all others except for the one of 



ref. j2l|, where we find that the difference is 2.5 times larger than the combined error. 



A. The RG-exponent yt^ 

We have determined the critical exponent yt^ from the slope of a renormalization group 
invariant quantity Ri at a fixed value i?2,/ of a second renormalization group invariant 
quantity R2 

S:=^ ^ cLy^^ , (29) 

where in our case Ri is either the Binder cumulant f/4 or the ratio of partition functions 
Za/Zp and we fix [Za/Zp) f = 0.3138. Fitting our data we have used the ansaetze 



^ = cL^'i , (30) 

S = cLy'^{l + dL-^) , (31) 

and 

5 = cL^*i(l + + eL-2) . (32) 

In table lylTl the results of fits for the slope of Za/Zp at Za/Zp = 0.3138 for the BCl model 
are given. We also performed joint fits of the BCl and the BC2 model. In table IVIIII we 
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TABLE VII. Fitting the slope of Za/Zp at Za/Zp = 0.3138 for the BCl model. 



Ansatz 




Vh 


d 


e 


xVdof 


m 


32 


0.71909(21) 






4.26/3 


m 


48 


0.71855(36) 






0.89/2 


m 


12 


0.71583(28) 


-0.148(6) 




20.91/5 


m 


16 


0.71700(38) 


-0.108(11) 




1.57/4 


m 


8 


0.71795(49) 


-0.026(21) 


-0.71(9) 


5.39/5 


m 


12 


0.71872(82) 


0.022(49) 


-1.01(30) 


3.75/4 



TABLE VIII. Fitting the slope of Za/Zp at Za/Zp = 0.3138 for the BCl and the BC2 model jointly. 



Ansatz 


Lrain 




di d2 


ei 


62 




m 


16 


0.71694(29) 


-0.109(9) 0.354(9) 






7.18/9 


m 


24 


0.71758(50) 


-0.084(21) 0.386(21) 






3.76/7 


m 


8 


0.71792(29) 


-0.027(11) 0.436(13) 


-0.71(4) 


-0.69(5) 


10.87/11 


m 


12 


0.71898(60) 


0.037(32) 0.505(32) 


-1.09(19) 


-1.12(19) 


4.86/9 



report results obtained from the slope of Za/Zp at Za/Zp = 0.3138. Also in the case of the 
slopes we observe that the amplitude of corrections is much larger for the BC2 than for 
the BCl model. Therefore using ansatz (130|) no acceptable fit can be obtained. 
Based on these fits we arrive at the final estimate 

yt, = 0.718(2) . (33) 

The central value and the error bar are chosen such that all the results of the fits reported 
in tables IVIII and IVIIII are covered. Our final estimate is also consistent with the results 
obtained from the slope of the Binder cumulant [4. We have also fitted our data for the 
Ising model with the ansaetze (130113 lll32p . The results obtained from the slope of Za/Zp 
for yt-^ are fully consistent with the estimate fl33|) . while those from the slope of f/4 are a 
bit larger. We conclude that leading bulk corrections oc L~'^ have only a small numerical 
effect in the case of the Ising model and can therefore be savely ignored in the case of the 
Blume-Capel model at Z) = 0.655, where leading corrections to scaling are suppressed at 
least by a factor of 30 compared with the Ising model. 



B. Surface susceptibilities and the RG-exponent y^^ 



Finally we have determined the exponent yh^ from the finite size scaling behavior of the 
surface susceptibilities Xii and Xi2- To this end we have computed 

X := x\R=Rf ^ L^'"-^-' (34) 

where we use (Za/Zp) f = 0.3138. We have fitted these quantities with the ansaetze 

X = cL^ (35) 
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TABLE IX. Fitting the surface susceptibilities xii and Xi2 at Za/Zp = 0.3138. Both the data of 
the BCl and the BC2 model are taken into account. Results for the parameters b and e are not 
reported to keep the table readable. 



Quantity Ansatz 


L'min 


X 


di 


d2 


xVdof 




ESI 


24 


1.2930(3) 


-0.09(1) 


1.13(1) 


10.37/7 




ESI 


32 


1.2932(4) 


-0.09(2) 


1.15(2) 


4.48/5 


Xii 


EZl 


16 


1.2935(4) 


-0.08(5) 


1.39(6) 


6.87/7 


Xii 


Ea 


8 


1.2933(6) 


-0.19(15) 


1.12(15) 


8.56/9 


Xl2 


ESI 


24 


1.2921(2) 


-0.36(1) 


0.02(1) 


4.01/7 


Xl2 


ESI 


32 


1.2922(3) 


-0.35(1) 


0.03(1) 


2.96/5 


X12 


E3 


16 


1.2932(4) 


-0.19(6) 


0.35(6) 


10.28/7 


Xl2 


Ell 


24 


1.2927(5) 


-0.25(12) 


0.16(13) 


3.10/5 


Xl2 


m 


8 


1.2935(2) 


-0.13(2) 


0.32(5) 


14.61/9 


Xl2 


m 


12 


1.2929(3) 


-0.30(4) 


0.02(9) 


7.83/7 


Xl2 


m 


16 


1.2928(6) 


-0.14(17) 


-0.06(25) 


3.19/5 



where x = 2yh^ — 2, 

X = cL-{l + dL-') , (36) 
X = cL"(l + dL^^) + h (37) 
where h is an analytic background and 

X = cL^(l + rfL-i + eL-2) + 6 . (38) 

It turns out that using the ansatz ( |35l) for none of the models studied an acceptable fit 
can be obtained. In table HXl results of joint fits of the BCl and the BC2 model using the 
ansaetze f p6|37|38p are reported. The variation of the results for the parameter x with the 
different ansaetze is of similar size as the statistical errors. The results obtained from xii 
are slightly larger than those from Xi2- Our final estimate x = 1.2929(10) covers all results 
reported in table IIXI Fitting the data obtained for the Ising model we get results for the 
parameter x that are essentially consistent with that obtained for the Blume-Capel models 
BCl and BC2. Therefore residual leading bulk corrections in the case of the Blume-Capel 
model at = 0.655 should not affect our final estimate of x. As our final estimate for the 
RG-exponent of the external surface field we quote 

2/;,, = 1.6465(6) . (39) 



IV. NUMERICAL RESULTS FOR THE SURFACE TRANSITION 

In the neighborhood of the special point, the transition line behaves as [l|-S] 

l^i.cl = c|t|* , (40) 

where $ = z/y^^ is the cross-over exponent. In the present section we shall compute the line 
of surface transitions numerically and check how well it is described by eq. ( HOj) . To this 
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TABLE X. Value of obtained by requiring that Za/Zp = 0.372884880824589... as a function of 
Lo and L. Throughout /? = 0.31. 



Lq L /?i,c 

10 8 0.6738189(96) 
10 10 0.6739017(78) 
10 12 0.6739022(65) 
10 20 0.6738941(40) 
10 100 0.6738924(95) 

5 100 0.6739308(83) 

6 100 0.6739055(86) 

7 100 0.6738980(85) 

8 100 0.6738761(85) 
12 100 0.6738845(82) 



end, we have only simulated the BCl model, i.e. the Blume-Capel model at D = 0.655 and 
Di = D2 = 0, since it has the smallest scaling corrections among the three models discussed 
in the preceding section. For given values of the bulk coupling /3 < /9c we have determined 
the critical surface coupling /3i c- To this end we have simulated films with an excess surface 
coupling /?! > at one surface and free boundary conditions, i.e. /32 = 0, at the other. 
For Lq ^ ^, where ^ is the correlation length of the bulk, the film should provide a good 
approximation of the semi-infinite system. We have determined Pi^c by requiring that the 
ratio of partition functions Za/Zp assumes the fixed point value of a square system with 
periodic boundary conditions in the universality class of the two-dimensional Ising model 



{Za/ZpY = 0.372884880824589... . (41) 



In the case of Za/Zp corrections to the fixed point value of the two-dimensional Ising univer- 



sality class vanish oc [23||. Note that in the case of f/4 and ^2nd/ L the analytic background 
of the magnetisation enters. Therefore, for these two quantities corrections vanish oc L~^/^. 

Similar to the previous section, we have updated the configurations by using a hybrid of 
the local heat-bath algorithm, the single-cluster algorithm and the wall-cluster algorithm. 
First we performed a series of simulations with varying values of Lq and L at /3 = 0.31, where 
^2nd = 1.04430(21). Our results are summarized in table |Xl We expect that corrections to 
the limit Lq,L — )■ 00 decay exponentially with increasing Lq and power-like as L increases. 
We see that for L = 20 and Lq = 10 the corrections are smaller than the statistical error of 
our numerical results. Assuming scaling, we chose for our simulations at larger values of P 
the lattice size such that Lq > 10^ and L = 2Lq. The numerical results of these simulations 
are summarized in table IXII In total these simulations took about 8 month of CPU time on 
a single core of a Quad-Core AMD Opteron(tm) Processor 2378 running at 2.4 GHz. 

In order to check the prediction (HUj) we have fitted our data with the ansatz 

/9i,c = /9i,sp + a U + ^Cif j (42) 
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TABLE XI. The surface coupling /3i^c at t-^e surface transition for various values of the bulk coupling 
p. In addition we give the second moment correlation length £,2nd of the hulk system, the linear 
lattice sizes Lq and L and the number of measurements (stat). 



2nd 



Lq L Stat/10^ 



/3 



l,c 



0.36 2.1272(4) 24 48 
0.378 4.1937(7) 50 100 
0.3822 6.0133(7) 70 140 
0.3842 7.9984(13) 100 200 
0.3859 12.1358(22) 150 300 
0.3865 15.6127(28) 180 360 
0.3869 20.0576(45) 250 500 



100 0.6394276(58) 
100 0.6093666(31) 
100 0.5967239(25) 
35 0.5883957(31) 
20.2 0.5785822(31) 
22.5 0.5738132(27) 
6.9 0.5698232(38) 



TABLE XII. Fitting the critical surface coupling /3i ^ with the ansatz H2\) . where m — 1 gives the 
number of correction terms that are included. For a discussion see the text. 



m 


1^1, min 




a 


C2 C3 




3 


0.36 


0.549350(9) 


0.51156(18) 


-11.27(7) 114.7(1.9) 


51.6/3 


3 


0.378 


0.549222(20) 


0.51490(50) 


-13.60(33) 266(21) 


0.77/2 


2 


0.378 


0.549454(7) 


0.50893(14) 


-9.45(4) 


154/3 


2 


0.3822 


0.549310(14) 


0.51236(32) 


-10.99(13) 


10.1/2 


2 


0.3842 


0.549240(26) 


0.51419(66) 


-12.15(39) 


0.00/1 


1 


0.3842 


0.549993(7) 


0.49483(12) 




877/2 


1 


0.3859 


0.549642(14) 


0.50214(28) 




28.8/1 


1 


0.3865 


0.549519(27) 


0.50492(59) 







where we have fixed $ = uyt^ = 0.63002 x 0.718. Our results are summarized in table IXlTl 
We see that without any correction, even the fit that includes only the three largest values of 
P that we have simulated has a very large x^/DOF and the result for /3i,sp is by far too large 
compared with that obtained in the previous section. Adding corrections terms, smaller 
values of /3 can be included into the fits. Also the value of /3i^sp gets closer to our result 
obtained above when adding further correction terms. However, the values of the correction 
amplitudes are large and rapidly increase with increasing order. We conclude that the be- 
havior of the surface transition seems to be consistent with the theoretical expectation (HOj) . 
However an accurate estimate of /3i,sp cannot be obtained from the numerical analysis of the 
surface transition. Finally we have performed a fit, where we have fixed Pi^sp as well as the 
crossover exponent Including all data with (3 > 0.378 we get 

/3i ,(/3) = 0.54914 + 0.51736(12) [t - 16.7(4)^^ + 752(101)t^ - 25015(6205)^^ ™''°'^^^ (43) 

where t = 0.387721735 - f3 with xV^OF = 1.23/2. This result might be helpful in future 
studies of the surface transition in the BCl model. 
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V. FILMS WITH (SP,+) BOUNDARY CONDITIONS 



Finally we have simulated the BCl model with special boundary conditions at one surface 
and + boundary conditions at the other, where + boundary conditions means that there 
are spins at xq = -^^o + 1 which are fixed to +1, which is equivalent to /i2 = /3. We have 
simulated at /3 = /3c and f3i = f3i^s = 0.549145, which was our preliminary estimate of 
Pi^sp at the time, when we started the simulations. We have measured the magnetisation 
profile. We have computed the Taylor expansion of the magnetisation profile in Pi around 
the simulation point Pi^s up to the second order. For each measurement we performed the 
following sequence of Monte Carlo updates: One sweep with the local heat-bath update, a 
cluster-update, one sweep with the local heat-bath update, and again a cluster-update. In 
these cluster-updates we fiip the sign of the frozen boundary and all spins of the cluster 
attached to it. Performing this cluster update twice, we are back to + boundary conditions 
at the surface. 

Mainly, we have simulated lattices with L = 4 x Lq. We have checked that at the 
level of our accuracy, deviations from the effectively two-dimensional thermodynamic limit 
L/Lo — 7- oo are negligible for this choice. To this end we have simulated L = 24, 32 and 
48 for the thickness Lq = 8. The results obtained for L = 32 and 48 are consitent among 
each other and those for L = 24 deviate only little. For L = 4 x Lq, we have studied films 
of the thicknesses Lq = 6, 8, 12, 16, 24, 32, 48, 64, 96 and 128. We have performed 10^ 
5 X 10^ 10^ 10^ 10^ 6 x lO^ 4 x lO^ 2.6 x lO^ 1.5 x 10^ and 7.7 x 10^ measurements for 
these lattice sizes, respectively. In total these simulations took a bit more than 2 years of 
CPU-time on a single core of a Quad-Core AMD Opteron(tm) Processor 2378 running at 
2.4 GHz. 



A. The magnetization at the boundary 



First let us discuss the behavior of the magnetisation mi at the boundary as a function 
of the thickness of the film. In terms of the reduced free energy per area it is given by 

dfiLo, t, h, ti, h) _ 1 1 E{s} exp(- + h J2x^ 



dhi 



L2 Z 



\x\,X2 I 



= mi 



(44) 



Note that the non-singular contribution to the free energy from the first surface does not 
feel the breaking of the symmetry by the second surface. Therefore it is an even function of 
hi and does not contribute to the partial derivative with respect to hi. Similar to eq. f lTSj) 
the singular part of the reduced free energy per area behaves as 



fs{Lo,t,h,ti,hi) = L^''+'g{...,hi[Lo/h,rn) 
Taking the partial derivative with respect to hi we get 



(45) 



mi 



dfs 



dhi 



dhi 



T-d+l 
^0 



t=h=ti=hi=0 



Vh-, 



cL 



t=h=ti=hi=0 
' 



(46) 
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TABLE XIII. Fitting the surface magnetisation mi at f3i = 0.54914. 



Ansatz 




X Ls d 




m 


8 


-0.35303(7) 0.908(3) 


31.87/6 


m 


12 


-0.35341(12) 0.938(8) 


14.37/5 


m 


16 


-0.35377(16) 0.975(14) 


3.13/4 


m 


6 


-0.35367(14) 0.988(14) 0.12(2) 


10.46/6 


m 


8 


-0.35396(21) 1.029(27) 0.19(5) 


6.97/5 



where denotes the partial derivative of g with respect to /ii[Lo///ii]^''i • 
We have fitted the magnetisation at the boundary with the ansaetze 

mi{Lo)=c{Lo + Lsr (47) 

with the free parameters c, Ls and x and 

mi(Lo) = c{Lo + LsT x (1 + d{L^ + Ls)'^) , (48) 

where we have included corrections oc Lq^ and oc i^o^, respectively. Using the Taylor- 
expansion, we have evaluated mi at /3i = 0.54914 and in order to estimate the effect of 
the uncertainty of our estimate of fii^sp also at /3i = 0.54916. Results of fits with the 
ansaetze (14 71148 p using the data for /3i = 0.54914 are given in table IXIlTl 

Repeating the fit with the ansatz ( H71) and -Lo,mm = 16 for the data taken at /3i = 0.54916 
we get X = -0.35280(16), = 0.935(13) and xVDOF = 2.98/4. Fitting the data for 
Pi = 0.54916 with the ansatz (gH]) and Lo,™„ = 8 we get x = -0.35288(18), = 0.965(22), 
d = 0.12(3) and x^/DOF = 6.08/5. We see that the error of the exponent x as well as that of 
Ls is dominated by the uncertainty of f3i^sp- As our final estimate we quote x = —0.3539(13) 
and hence 

Vh, = 1.6461(13) (49) 

which is fully consistent with but less precise than the estimate obtained above, eq. (!39ll . 
For the effective shift in the thickness of the lattice we quote Ls = 1.0(1). Note that Ls 
should be equal to the sum of the extrapolation lengths at the two surfaces. For + boundary 



conditions we find in ref. j2J] the result Ig^ = 0.96(2) for the + boundary conditions. Hence 
hx,sp ~ for the BCl model. This is consistent with our observation above that for the BCl 
model corrections oc Lq ^ have a small amplitude. 

We also have performed fits with the ansaetze (I47ll48p . where f3i^sp is a free parameter. 
To this end we have used the Taylor series of the surface magnetisation mi in /3i around the 
simulation point. Such fits suggest Pi^sp = 0.54916(4), where the error is larger than that of 
our final estimate fl27j) . Correspondingly, also the estimate of i/hi cannot be improved this 
way. 

B. The universal scaling function of the magnetisation profile 

The magnetisation profile at the critical point behaves as [l|-0] 

m{z) = c ij{z/Lo) (50) 
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0.2 0.4 0.6 0.8 1 

z/Lo,eff 



FIG. 2. We plot m{z)Lj^j J as a function of z/Lg^e//) where z is the distance from the boundary 
with special boundary conditions. The data are taken at the bulk critical point. At the first 
boundary special boundary conditions and at the second + boundary conditions are applied. We 
give data for the thicknesses Lq = 8, 32 and 128. The line given for Lq = 128 is obtained from 
a linear interpolation of the data points. We abstained from giving data for more thicknesses to 
keep the figure readable. Notice that the statistical errors are smaller than the symbols and the 
line shown. For a detailed discussion see the text. Colour online. 

where z gives the distance from the boundary. The universal function ip depends on the 
surface universality classes of the boundary conditions at the two surfaces of the film. In 
the neighborhood of the surface, the magnetisation follows a power law 

m{z) = az-^/'^+'^-^-y'^^ . (51) 

Note that the exponent —P/v + ci — 1 — = —0.1646(7) is negative and therefore the 
magnetisation is actually decreasing with increasing distance from the surface. This is in 
contrast to ordinary boundary conditions, where i/hi = 0.7249(6) and therefore the 
magnetisation is increasing with the distance. Note that from scaling relations it follows 
that /3/z/ = (1 + ri)/2, where 1] = 0.03627(10) for the three-dimensional Ising universality 
class [5|]. 

Using the Taylor-expansion up to second order we have computed the magnetisation 
profile for /3i = 0.54914, which is our estimate of /3i,sp. In figure[2]we have plotted m{z)L^^^^j 
as a function of z/LQ^f.ff- In order to take corrections oc Lq ^ into account, we have replaced 
in eq. ( I5U]) the thickness Lq by the effective one i^o,e// = lex,sp + lex,+- The distance from the 
boundary is given hy z = XQ — l/2 + lex,sp- As discussed above, the extrapolation lengths take 
the values lex,+ = 0.96 and lex,sp = 0. We see that the data for the three different thicknesses 
fall nicely on top of each other. The same holds for thicknesses not plotted here. Therefore 
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the curve obtained from Lq = 128 should be a good approximation of the universal scaling 
function iIj{z/Lo). We observe that, as discussed above, for small z the magnetisation is 
indeed decreasing with increasing z. At z / Lq ^ 0.31 a very shallow minimum is reached. 



VI. SUMMARY AND CONCLUSIONS 



In this work we have studied the special point in the phase diagram of a semi-infinite 
system. We have accurately estimated the surface RG-exponents yt^ and that govern, 
along with the bulk RG-exponents yt and yt, the behavior of surface quantities at the 
special point. To this end, we have simulated the improved Blume-Capel model with two 
different surface interactions and the spin-1/2 Ising model on the simple cubic lattice. We 
have determined the surface RG-exponents from a finite size scaling analysis at the special 
point. In the presence of surfaces typically corrections oc Lq ^ appear, where Lq is the linear 
extension of the finite system. Fitting data it is difficult to disentangle these from the leading 
bulk corrections oc Lq'^, where u = 0.832(6) [sj. Therefore it is helpful that in the improved 
Blume-Capel model the amplitude of the leading correction to scaling vanishes. 

In table IXTVl we compare our results with those obtained in previous Monte Carlo (MC) 
studies and computed by field theoretic methods. In previous MC studies only the spin-1/2 



Ising model on the simple cubic lattice had been simulated. Only in ref. [21| results for 
the RG-exponents are quoted. However all authors give results for the exponents /3i and 
most for $. In table IXTVl we have converted these estimates by using the scaling relations 
f3i = u{2 — and $ = uyt^, where we have used z/ = 0.63002. In the case of yh^ we 



observe that the estimates of refs. |2ll . |25| deviate by 9 and 6.6 times the combined error 



from our result. Notice that such deviations can be caused by corrections that are not taken 
into account in the ansatz. For a general discussion of this problem see section Hi CI The fits 



performed in ref. [2l| are in particular prone to this problem, since in the ansatz corrections 



oc Lq^ are included while e.g. corrections oc L^"^^ are missing. In the present work, we 
tried to estimate such systematical errors by comparing the results obtained by fitting with 
ansaetze of a varying number of correction terms. In the case of yt^ we see a quit e g ood 
agreement with the results of refs. (2^, 2l||. In contrast, the estimates of refs. [l8[ 19, 27 



are clearly larger than ours. 

The surface critical exponents for the special point have been computed up to second 



order in the e-expansion. Direct results have been obtained for the exponents t^h [28| and 
$ j29|. Starting from these two exponents and the bulk exponents, the remaining ones can 
be computed by using scaling relations. Naively inserting e = 1 one gets rj\\ = —0.30 and 
$ = 0.68 for three dimensions; see table VI of ref. [2|. In table IXTVl we give yti = {d — r]\\)/2 
and yt^ = where we have used z/ = 0.63002. The result for yh^ is in quite good 

agreement with those from Monte Carlo simulations. In contrast, the estimate for yt^ is by 
far too large. The authors of ref. [30| have pointed out that by resummation of the series. 



e.g. by Fade approximants this discrepancy can be resolved. The authors of ref. [30| have 
computed the surface exponents in a two-loop calculation in three-dimensions fixed. The 
numerical estimates quoted by the authors are obtained by using Fade approximants. The 
numbers given in table IXIVI are again obtained by using yh^ = {d — rj\\)/2 and yt^ = $/z/, 
where v = 0.63002. The results for yh^ and in particular for yt-^ deviate clearly from those 
obtained from Monte Carlo simulations. 

In addition to the estimates of the RG-exponents, our finite size analysis provides us with 
accurate estimates of the surface coupling f^i^gp at the special point for the three models that 
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TABLE XIV. Results for the surface RG-exponents of the special phase transition in the three- 
dimensional Ising bulk universality class. 



Ref. 


year 


Method 




Vh 


[18j 


1984 


MC 


1.72(4) 


0.89(6) 


[19] 


1990 


MC 


1.71(3) 


0.94(6) 


[22 


1992 


MC 


1.65 


1.17 


[20] 


1992 


MC 


1.624(8) 


0.732(24) 


[25] 


1995 


MC 


1.623(2) 




[26] 


1998 


MC 


1.635(16) 




[2]J 


2005 


MC 


1.636(1) 


0.715(1) 


this work 


2011 


MC 


1.6465(6) 


0.718(2) 


[2,28, 29] 


1981 


e-exp, naive 


1.65 


1.08 


[29, 30] 


1998 


e-exp, resummed 




0.752 


[30] 


1994 


3D-exp, resummed 1.583 


0.856 



we have simulated. These estimates can be used in future studies of thin films. In particular, 
we intend to study the thermodynamic Casimir force in thin films in the neighborhood of 
the special point. 

Furthermore we have studied the behavior of the surface transition in the neighborhood 
of the special point. Our numerical results follow the theoretical expectations. However the 
special point cannot be located accurately with such an approach. 

Finally we have simulated a film with symmetry breaking boundary conditions at one 
surface and special boundary conditions at the other. The behavior of the magnetisation 
at the surface with special boundary conditions follows a power law. Its exponent can be 
expressed in terms of the RG-exponent yh^- The analysis of the data fully confirms our 
estimate of yh^ given in table IXIVI The magnetisation profile follows a universal function. 
This theoretical expectation is fully confirmed by the nice collapse of data that we observe 
for a large range of thicknesses Lq of the film. An interesting feature of the magnetisation 
profile is that at the surface with special boundary conditions, which do not break the 
symmetry, the magnetisation is decreasing with increasing distance from the surface. 
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